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ABSTRACT 

Quasar optical monitoring campaigns are reaching a new standard of quality as a result of long 
term, accurate CCD observations. In this paper we review the basic Poissonian formulation of quasar 
variability, using it as a mathematical tool to extract relevant parameters such as the energy, rate and 
lifetimes of the flares through the analysis of observed light curves. It is shown that in this very general 
framework the well established anti-correlation between variability amplitude and wavelength can only 
be understood as an effect of an underlying spectral component which remains stable on long time-scales, 
and is redder than the variable component. The formalism is applied to the B and R light curves of 42 
PC quasars collected by the Wise Observatory group (Giveon etal. 1999). Variability indices for these 
data are obtained with a Structure Function analysis. The mean number of living flares is constrained 
to be in the range between A/" ~ 5 and 100, while their rates are found to be of order v ~ 1-100 
yr~^. Monochromatic optical flare energies E\ ^ ^.O'*^"'*^ ergA~^ and life-times r of 0.5 to 3 yr are 
derived. Lower limits of typically 25% are established for the contribution of a non-variable component 
in the R band. The substantial diversity in these properties among quasars invalidates simple versions 
of the Poissonian model in which flare energies, lifetimes and the background contribution are treated as 
universal invariants. Light Curve simulations confirm the applicability of the method. Few significant 
correlations between variability indices and multi- wavelength properties of the quasars exist, confirming 
the results of Giveon et al. The good correlation between the EW(H/3) and the long term variability 
amplitude is interpreted in a scenario where only the variable component participates in the ionization 
of the line emitting gas. This idea is consistent with the observed trends of the variability amplitude 
with A, EW(HcII) and the X-ray to optical spectral index. The parameter estimates derived under the 
framework of Poissonian models are applicable to several scenarios for the nature of quasar variability, 
and can help guiding, testing and discriminating between detailed physical models. 

Subject headings: quasars: general - galaxies: active - galaxies: Seyfert - galaxies: nuclei - methods: 
statistical - methods: analytical 



1. INTRODUCTION 

Variability is one of the defining properties of Active 
Galactic Nuclei (AGN) and a potentially powerful dis- 
criminant among different scenarios for the physics of the 
central engine. Optical-UV variability studies in the past 
decade have put a lot of effort into investigating the effects 
of continuum variations upon the emission lines, leading to 
substantial progress in the diagnostics of the properties of 
the Broad Line Region (Netzer & Peterson 1997). Yet, no 
comparable progress has been achieved in understanding 
the origin of the continuum fluctuations. 

In this paper we discuss quasar variability in the context 
of Poissonian models. By Poissonian we mean any physi- 
cal system in which the variations are due to the stochastic 
superposition of independent flares, occurring at a given 
mean rate but randomly distributed in time, a process also 
called "christmas-tree" , "shot-noise" , "discrete events" or 
"subunits" model. This approach has two key advantages: 



(1) It provides a simple mathematical framework which 
can be applied to observed light curves to constrain the 
most relevant flare properties, such as energy, time-scale 
and rate (e.g., Cid Fernandes, Aretxaga & Terlevich 1996). 

(2) It encompasses a large class of physical models for 
AGN. For instance, scenarios as diverse as accretion-disk 
instabilities of several kinds (Haardt, Maraschi & Ghis- 
ellini 1994; Kawaguchi etal. 1998), disruption of stars in 
the gravitational field of super-massive black-hole (Peter- 
son & Ferland 1986; Ayal, Livio & Piran 2000), stellar 
collisions (Keenan 1978; Courvoisier, Paltani & Walter 
1996), supernovae (Aretxaga & Terlevich 1994; Aretxaga, 
Cid Fernandes & Terlevich 1997), and even extrinsic vari- 
ability models such as micro-lensing (Hawkings 2000) all 
share a common Poissonian nature. It is the combina- 
tion of these two points that makes the Poissonian inter- 
pretation of AGN variability attractive, as its application 
to observed data may provide feasibility tests on different 
theories. 
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The early realization that most luminous sources tend 
to vary less (Uomoto, Wills & Wils 1976; Pica & Smith 
1983) qualitatively favored a Poissonian interpretation. 
However, these and subsequent quasar monitoring stud- 
ies found that the slope of the fractional variability 

{6 = a{L)/L) versus mean luminosity relation, 5 (x L , is 
shallower than the (3 = —1/2 value predicted in a sim- 
ple Poissonian model (Cristiani, Vio & Andreani 1990; 
Trevese et al. 1994; Hook et al. 1994, Cristiani et al. 1996, 
Paltani & Courvoisier 1997), which lead some of these 
studies to rule out such models. Some studies, however, 
find slopes consistent with (3 = —1/2 (Cid Fernandes et al. 
1996; Garcia etal. 1999), while others even question the 
very existence of a correlation (Bonoli etal. 1979; Net- 
zer & Sheffer 1983; Giallongo, Trevese & Vagnetti 1991; 
Lloyd 1984; Cimatti, Zamorani & Marano 1993; Netzer 
etal. 1996). 

The slope of the variability-luminosity relation has thus 
been a controversial issue, with conflicting results reported 
in the literature. Selection effects and different cover- 
ages of the luminosity-rcdshift plane are likely causes for 
these discrepancies, as discussed in Hook etal. (1994). 
Garcia etal. (1999) argued that, although the observed 
variability-luminosity relation in their sample is shallower 
than (3 = —1/2, good agreement between data and model 
is achieved after taking into account the increase in vari- 
ability with frequency. Moreover, Cid Fernandes etal. 
(1996) argued that the photometric (at least for pho- 
tographic data) and sampling uncertainties, along with 
wavelength/rcdshift effects, propagate to a poorly defined 
variability-luminosity relation. Furthermore, they showed 
that a slope of -1/2 is only expected in the simplest of 
Poissonian models, in which the flare energy, time-scale 
and background contribution are held fixed as universal 
constants for all objects. Discarding Poissonian models on 
the basis of /3 7^ —1/2 is thus both risky and an oversim- 
plified interpretation of Poissonian processes. 

Quasar variability studies are reaching a new level of 
quality due to the efforts of several groups which engaged 
into long term, differential CCD photometric and spectro- 
scopic monitoring programs (Borgeest & Schramm 1994; 
Netzer et al. 1996; Sirola et al. 1998; Giveon et al. 1999— 
hereinafter G99; Garcia etal. 1999; Kaspi etal. 2000). 
These studies represent an enormous improvement over 
pre-existing variability databases, both in sampling and 
photometric quality. While most previous studies (by ne- 
cessity) based their analysis on properties of the ensemble 
of objects, these new data allow the study of quasar vari- 
ability properties on an object-by-object basis. With the 
continuation and steady improvement of these programs, a 
clearer picture of the phenomenology of quasar variability 
will eventually emerge. 

In this paper we review the formal relations between 
observable variability properties and the parameters in a 
general Poissonian model (§2). This approach is here seen 
as a valid step towards a physical understanding of the 
nature of AGN variability, which is particularly relevant 
given the current lack of an accepted paradigm for this 
ubiquitous phenomenon. Emphasis is put on the use of 

^To be precise, if p(x) is the probability density of x, one finds 



multi-wavelength data to constrain the basic model pa- 
rameters. In §3 we apply the formalism to the photomet- 
ric (B and R) quasar monitoring data collected over the 
past decade in the Wise Observatory (G99), by means of 
basic light curve statistics and an analysis of the individ- 
ual Structure Functions. Simulations are used to verify 
the consistency of the results. A discussion on the use of 
higher moments of the data (skewness and kurtosis) is also 
presented (Appendix A). In Section 4 we discuss our re- 
sults and presc!nt a correlation analysis of the variability 
with other properties, as well as possible interpretations in 
the context of AGN models. Finally, in Section 5 we sum- 
marize our main results and outline prospects of future 
work. 

2. FORMALISM 

Poissonian models for AGN variability have been stud- 
ied both qualitatively and quantitatively, and for wave- 
lengths across the electromagnetic spectrum, from radio 
(e.g.. Dent 1972), to optical-UV (e.g., Cid Fernandes etal. 
1996; Garcia etal. 1999) and X-rays (e.g., Lehto 1989; Al- 
maini etal. 2000). In this section we follow the detailed 
formulation of this problem developed by Cid Fernandes 
(1995). The resulting formulae rest upon basic probability 
and time-series theory (e.g., Papoulis 1965). 

In a Poissonian scenario the instantaneous monochro- 
matic luminosity L\{t) is due to the superposition of a 
variable component, V\{t), and an underlying background 
component C\: 



Lx{t) = Vx{t) + Cx 



(1) 



where Vx{t) is made by the superpositions of flares lx{t) 
with random "birth-dates" tf. 



Vx{t) = ^lxit-U) 



(2) 



The first two moments (mean luminosity and relative 
variability) of the variable component are 



Vx = uxEx 



<j{Vx) _ 1 
Vx {yxTxY'^ 



(3) 



(4) 



In these expressions vx is the mean rate of flares. Ex is 
the monochromatic energy of individual flares, and tx is 
the flare "life-time" , defined by 



TX 



El 



film 



(5) 



rx 



It can be shown that these relations also hold when one 
allows for the (likely) possibility that the flares within a 
given object are not all identical, i.e., when lx{t) = lx{t, x), 
where x denotes a general set of parameters (size, density, 
cooling time, ...). The only modifications in this more 
general case is that Ex above means the average Ex{x) 
over the probability distribution of x, and similarly for 



[(t, x)p(x)dxdt 
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The total mean luminosity and relative variability are 
both affected by the underlying constant component, turn- 
ing eqs. (3) and (4) into 



Lx = y\Ex + Ca 



(6) 



(7) 



where v\ = V\/L\ is the fraction of the mean luminos- 
ity which is actually due to the variable component. Al- 
lowance for the possible contribution of a non-variable 
component is essential, though it has often been forgotten 
when discussing Poissonian models. Physically, C\ can 
be associated with several sources. These can be either 
around the nucleus and extrinsic to the variability gen- 
eration process, as the host galaxy stellar contribution, 
or, more interestingly, a part of the variable continuum 
source which remains stable over long time-scales, such as 
the non-flaring part of an accretion disk. 

The Poissonian model thus involves four basic parame- 
ters: the rate [vx), energy {Ex) and life-time (ta) of the 
flares, plus Cx- The shape of the flares may be regarded 
as a further degree of freedom, but it has little effect upon 
the analysis presented here. Equations (6) and (7) relate 
these parameters to just two observables. Even estimat- 
ing Ta through a Structure Function (SF) analysis, one 
has a non-closed system, with three observables and four 
variables. Higher moments of the light curve could in prin- 
ciple be used as further constraints (see Appendix A), but 
these are so badly affected by sampling uncertainties that 
presently they do not provide useful constraints. Another 
piece of information that can help constraining the model 
parameters is that the constant component has to be at 
most as strong as the observed minimum in the light curve: 



< Ca < Lx, 



(8) 



It is reasonable to presume that the same flares are seen 
across a narrow spectral band, such as the optical UV. 
Empirical support for this hypothesis comes from the high 
similarity between the continuum light curves in different 
bands (e.g., Cutrietal. 1985; Kroliketal. 1991; G99). In- 
deed, the whole method of reverberation mapping relies 
on an equivalent hypothesis, namely, that the fluctuations 
in the ionizing continuum can be mapped by those of the 
optical-UV continuum (Peterson 1993). We therefore may 
write z^A = and tx = t. The developments below could 
also be made allowing for A dependent time-scales, but we 
shall adopt the simpler constant ta scenario. 

A first consequence of this assumption is that one can 
isolate the spectral shape of flare energy directly from the 
spectral behavior of the standard deviation <jx'- 



Ex 



1/2 



ax oc (TA 



(9) 



It is well established that the amplitude of the vari- 
ations increase towards shorter wavelengths (Cutri et al. 
1985; Edelson, Krolik & Pike 1990; Kinney etal. 1991; 
Cristiani etal. 1997; Di Clemente etal. 1996), with most 
sources becoming bluer as they brighten, which immedi- 
ately tells us that the variable component is blue. This 



fact also has the very interesting consequence (see eq. [7]) 
that vx must vary with wavelength, which can only he un- 
derstood invoking an underlying component. 

Another way of seeing this is to rewrite (7) as 5x = 
VxM~'^/'^, where Af = vt is the mean number of living 
flares at any time. Since we are assuming that the same 
flares are seen in different wavebands, the fact that 5x de- 
creases towards the red can only be made consistent with 
a Poissonian scenario if the relative contribution of the 
background, 



CA 



^ = 1 
Lx 



increases with A. Of course, the decomposition of optical- 
UV spectra of AGN into a variable plus a constant compo- 
nent has been proposed before, both on observational and 
on theoretical grounds. Here we proved that the spectral 
behavior of Sx implies the existence of a constant source if 
one is to keep within the framework of Poissonian models. 

One can construct families of possible spectra for the 
constant component by writing 



Cx 



O-A 



(10) 



The simultaneous analysis of the light curve statistics 
in different wavebands thus has interesting consequences, 
but we still do not have a closed system. Except for the 
life-time t, which can be estimated through a structure 
function analysis, there is no way to determine absolute 
values for i/, Ex or Cx- Nonetheless, the wavelength infor- 
mation, coupled with the condition that Cx must satisfy 
(8), can yield improved constraints upon Af. One may 
combine the above relations to obtain 



MA 



5x 



<J\f< 



Sx 



(11) 



where /iA = Lx^min/ Lx- Both lower and upper limits can 
be made more stringent by considering the whole wave- 
length information available. We can thus define Mmin by 
using the maximum value of the lower limit above, and 
conversely for TVmax- This range of allowed A/" translates 
into corresponding ranges for v, Ex and Cx- Exactly how 
stringent A^min is depends on how close LA,min gets to Cx ■ 
As the chances of the light curve reaching the background 
level decrease with the increasing superposition of events, 
one expects the A^min limit to get progressively less strin- 
gent as M increases. The quadratic dependence on the 
observed quantities also conspires to broaden the range of 
M. 

This very general and straightforward formalism can be 
applied to several variability data sets. In the next section, 
we apply it to one of the best sets of quasar optical light 
curves presently available. 



3. ANALYSIS OF THE WISE OBSERVATORY QUASAR 
LIGHT CURVES 

3.1. Description of the data 
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Givcon ct al. (1999) have presented the resuhs of a long 
term B and R photometric monitoring of 42 opticaUy se- 
lected nearby quasars from the Palomar Green (PG) sam- 
ple. The observations were collected with the Wise Obser- 
vatory Im telescope over the 1991 to 1998 period. Spcc- 
trophotometric measurements at Wise and Steward Obser- 
vatories were used to complement the light curves of 13 of 
the objects. Previous results on this and related monitor- 
ing campaigns have been pubhshed in Maoz etal. (1994), 
Netzer etal. (1996), and a spectroscopic study of a sub- 
sample of the objects has been recently concluded (Kaspi 
etal. 2000). The reader is referred to G99 for details of 
the observations. The following quantities are medians 
over the sample: Uots = 32 observations per object, rest- 
frame sampling interval of 33 days, rest-frame light curve 
span of 5.9 years, photometric uncertainty = 0.015 mag 
(in B), z = 0.16; Mb = -23.3. 

The apparent magnitude light curves were converted 
to monochromatic luminosities Lx at 4400 (B) and 6400 
A (R) using the same cosmological parameters as G99 
{Ho = 70 kms-iMpc-i, qa = 0.2, A = 0). The K- 
corrcction was also performed as in G99, assuming a 
power-law spectrum, whose slope is defined by the me- 
dian B — R color. This was a minor correction because 
of the low redshifts in this sample. Galactic extinction 
corrections were made with the values oi As extracted 
from NED'', and the extinction law of Cardelli, Clayton & 
Mathis (1989, with Ry — 3.1), but were mostly negligible 
{Ab < 0.17, median = 0.03) because of the high galactic 
latitude of the objects (Schmidt & Green 1983). Emission 
lines contribute ^ 10% to the B and R band huninosities 
of PG quasars (G99); their effect upon the present analysis 
is also negligible. 

The analysis presented below was carried out with the 
resulting LB{t) and Lji{t) light curves. We note that the 
estimates of A/", ly, r and the background fractional con- 
tribution are cosmology, extinction and K-correction inde- 
pendent. Only the absolute values of E\ and C\ depend 
on such factors. 

3.2. Structure Function Analysis: Estimating time-scales 

and the asym,ptotic variance 

The most commonly employed tool to extract informa- 
tion on the variability time-scales out of quasar light curves 
is the Structure Function (Simonetti, Cordes & Heeschen 
1985; Hook etal. 1994; Cristiani etal. 1996). The SF is 
defined by the mean of [L{t + At) — L{t)Y over a light 
curve. For a Poissonian sequence of flares the SF is known 
to be simply proportional to the SF of a single isolated 
flare (Appendix B), and therefore only has structure for 
time-scales shorter than the flare duration. For large At 
the SF converges to twice the intrinsic variance of the pro- 
cess: SFx{M) 2al = 2EIv/t (eq. [7]). This happens 
because variations on such long time-scales effectively cor- 
respond to independent samples of the L\{t) process, as 
the light curve loses memory of its past. 

The SF provides estimates for both the variability time 
scale and amplitude. We use the notation a\ gp to dis- 
tinguish the asymptotic variance derived from the SF 
from that computed directly from the light curve a\ = 

^The NASA/IPAC Extragalactic Database (NED) is operated by 
contract with the National Aeronautics and Space Administration. 



L\ — Lx — e\ (where the last term corrects for the small 
effects of photometric error), which may be somewhat un- 
derestimated for light curves not much longer than the 
correlation time-scale since they do not sample the whole 
power contained in the fluctuation power spectrum. Like- 
wise, 5x,SF = o'x^sf/Lx denotes the long term net vari- 
ability. 

The SFs were fit with theoretical functions correspond- 
ing to different flare evolution models. In all flts we are 
interested in just two quantities: The flare life-time r, de- 
fined by (5), and gp. Four models for l{t) were explored: 
(1) square; (2) exponential; (3) symmetric triangle and 
(4) asymmetric triangle. Equations and plots for the cor- 
responding SFs are displayed in Appendix B. These are 
clearly toy-models for the radiative evolution of physical 
flares, but as shown in Appendix B (Fig. BIO) the SF is 
not very sensitive to l(t). While it is unfortunate that lit- 
tle information about l{t) can be retrieved with this tech- 
nique, this does make the estimates of r and 

c^x SF ^ore 

robust. This is confirmed by the fits, which showed that 

all l(t) models produce very similar parameters (with the 
exception of exponential flares, which in some cases do 
not converge and are generally poorer). For this reason, 
we present only the results for square flares. 

Let Li [i = 1, ... TV) be the luminosity of a quasar at 
an epoch ti in a given spectral band. The SF of this light 
curve may be defined as follows. If ti and tj are two dis- 
tinct epochs of the light curve, and ti > tj, then 

= s{At,,) = [L, - L, f - (e, + Cj f (12) 

where At ^ = {ti—tj)/{l-\-z) is the interval between these 
two epochs in the quasar rest-frame. In this expression, 
the observational errors e are subtracted in quadrature 
from the luminosity difference in order to "remove" the 
observational noise from the SF. The data to be fitted, 
hence, are the set of pairs (Aty , Sjj). The SF parameters 
were estimated by direct minimization of the sum of the 
square residuals between data and model. Note that no 
binning was applied to data before the fitting. All but 
the fourth model in Appendix B have two parameters: a 
time scale r and the asymptotic limit SF{oo) = 2<7\ gp 
(the fourth model in Appendix B has two time-scales). 
The errors in the parameters were obtained by bootstrap. 
In this technique an observed distribution with Uobs data 
points is resampled many times, each time picking a ran- 
dom set of Uobs data points (allowing for repetition) out 
of the original ones. The underlying hypothesis is that 
the actually observed values trace the distribution of the 
measured quantities. For each resampled data set the pa- 
rameters are fitted, and the parameter errors are estimated 
from the variance of the values estimated in all resampled 
data sets. Initially the B and R SFs were fitted separately. 
A global estimate of r combining the B and R data was 
then performed, fitting both SFs with the same value of r 
but different a% gp and a\ gp . 

B and R SFs for five objects are illustrated in Fig. 1 
along with their individual and global fits. The SFs in 
this figure were computed binning the squared luminosity 
differences in At such that each bin contained 25 points. 

Jet Propulsion Laboratory, California Institute of Technology, under 
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Fig. 1. — Illustration of the Structure Functions and their fits for 5 quasars. Left and right panels show the results for B and R bands 
respectively. Dotted lines correspond to the individual SF fits, while solid lines indicate the SF resulting from the combined B and R fits. 
The {L{t + At) — L{t))^ differences were grouped in bins of 25 points and the error bars in the SF were computed with the bootstrap method. 
This was done for plotting purposes only, since the fits do not involve binning the SF. All fits shown correspond to square flares. The SFs are 
normalized to twice the observed light curve variance (2(t^) for convenience. 



This was done just for plotting purposes, as the fitting pro- 
cedure does not involve binning. The error bars in Fig. 1 
were computed bootstrapping the light curves 1000 times. 

Though the fits were satisfactory for most objects (e.g., 
PG 0838 and 1307 in Fig. 1), very often the SFs exhibit 
complex shapes with ups and downs (e.g., PG 1613 and 
2251). Such a non monotonic behavior does not rule out 
the hypothesis that the light curve is made up by the Pois- 
sonian superposition of flares with simple profiles, as large 
departures from the predicted SF also occur for simulated 
light curves when these are sampled a finite number of 
times, as will be explained in Fig. 4 and §3.4.1. Theoreti- 
cal SFs are derived under the assumption of infinitely long 



light curves, and it would be naive to expect observed SFs 
to exhibit the smooth, well behaved shapes predicted by 
theory. This point has to be kept in mind when evaluating 
the fits performed in this section. One must therefore exert 
caution when interpreting the fit parameters, particularly 

T. 

Results of the fits are presented for all 42 sources in 
Table 1, along with other useful light curve statistics. 
Though in many cases the life-times inferred from both 
wavelengths agreed well, we find a large scatter around 
the tb = tr fine (see also G99). For the reasons discussed 
above, we believe this is more likely an artifact of the fi- 
nite span of the light curves than a real effect, and in what 
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follows we shall only refer to the values of t obtained with 
the combined B and R fits. The asymptotic net standard 
deviations Sb,sf and Sr^sf in Table 1 also correspond to 
those obtained in the combined fits, but we remark that 
these were very similar for the individual and global fits. 

For four objects, PG 0026, 1100, 1427 and 1626, we find 
that the SF does not converge over the time span of the 
observations. All these quasars present rising or decreas- 
ing trends over all the length of the observations (Fig. 2 
of G99). This suggests cither long (> 5 yr) time scales, 
an underlying slowly varying component, or else may be 
indicative of a non-statistically stationary process. In ad- 
dition, for PG 1114 and 1512 we find time-scales close to 
the length of the light curve. These six objects were dis- 
carded from further analysis. 

Fig. 2 shows how the net variability measured directly 
from the light curve compares with the SF-determined 
long term net variability. One finds that Sx,sf is some 
17% larger than S\ for both B and R, confirming that 
the latter slightly underestimates the long term variabil- 
ity. Also shown in Fig. 2c is a comparison of the net vari- 
abilities in the B and R bands, to illustrate the fact that 
Sb,sf > Sr^sf by typically ^ 30%. 

The diversity of time-scales obtained in this analysis 
(from 0.5 to 3 yr) is consistent with the results of G99, 
who found a wide range of zero-crossing auto-correlation 
time-scales. In fact, our values of r are similar to theirs, 
but with a substantial scatter. In the context of Poissonian 
models, this diversity reinforces the suspicion (§1) that 

— 1/2 

quasar flares are not all the same, and that the 6 (x L 
may not be universal. Indeed, in this more general sce- 
nario, the "generalized Poissonian model" of Cid Fernan- 
des etal. (1996), every quasar lies in one of a family of 

1/2 

6 (X L laws (see Fig. 6). Members of each such family 
have difi^erent rates, but same energy, life-time and back- 
ground contribution (eqs. [6] and [7]). Also, this intrinsic 
diversity casts doubts upon the meaning and usefulness of 
ensemble SFs. 

3.3. Constraints on the energies, rate and background 

contribution 

The asymptotic net variabilities dB,SF and Sr^sf ob- 
tained above, along with the values of L\ and fix (listed 
in Table 1), can be directly applied on eq. (11) to constrain 
the allowed range of J\f, cb and cr. We remark that this 
is nearly independent of the SF analysis, since S\^sf is 
not too difi^erent from Sx (Fig. 2) and r is not employed 
in these constraints. The results of this calculation are 
presented in Table 2. 

Except for 3 out of 36 quasars, the upper limit Afmax 
is always given by the B band, since, as already noted by 
G99, the variations there are larger than in R (Fig. 2c), in 
agreement with the general tendency of AGN. As antici- 
pated (§2), the ranges of the obtained TVs are wide (factor 
of ~ 20) . Still, they provide useful estimates of a physically 
meaningful quantity in Poissonian models. From Table 2 
one concludes that the typical number of living events in 
the G99 quasars is of order 5-100. 

A corollary of Mmax being defined by the B light curves 
is that non trivial lower limits for the background compo- 
nent can only be obtained for the R band for the major- 



ity of objects (Table 2). In some cases, like for PG 1309 
and 1354, one finds lower limits of as much as 40% in R, 
indicating a substantial contribution from an underlying 
constant spectral component. The upper limits for cr are 
given by jiR, tabulated in Table 1, which for these two 
quasars are 90 and 86% respectively. 

3.3.1. Energies and rates 

The N and cx limits above make no explicit use of the 
life-times. Estimates for the allowed range of flare rates 
and energies require knowledge of r. These arc listed in 
Table 2 for the life-times found in §3.2. Rates between 
1 and a few hundred yr^^, and monochromatic energies 
from 3 X 10^^ to 10''° ergA"! are found for Eb and Er. 
Note that the energies for individual quasars are more con- 
strained (typically to within a factor of 4) than either v or 
AT, since they only depend on the square root of v. Fur- 
thermore, inspection of Table 2 shows that there is not 
a single value of either Er or Er that is simultaneously 
compatible with all lower and upper limits for the G99 
quasars. This, again, points to a diversity of flare proper- 
ties in quasars. 

3.4. Light curve simulations 

The estimates of the basic Poissonian parameters pre- 
sented here rely on just the mean and minimum luminosi- 
ties plus the SF-based estimates of the flare life-time and 
asymptotic variance, this latter quantity being little differ- 
ent from the light curve variance. The actual light curves 
contain much more information, as they are defined by a 
particular realization of birth-dates of the flares and their 
detailed radiative evolution. Retrieving this information 
from these observations is however an impossible task. We 
have experimented using higher moments (skewness and 
kurtosis) of the light curves as further constraints, but 
this proved fruitless in practice (Appendix A). 

Since our estimates are based on so little information, 
there is no guarantee that a superposition of flares with 
the inferred properties would bear any morphological re- 
semblance to the observed light curves. In order to ver- 
ify whether Poissonian models within the bounds defined 
in §3.3 can produce quasar-looking light curves we have 
performed a series of Monte Carlo light curve simulations. 
Our goal here is to broadly assess the "morphological com- 
patibility" of models and data in a qualitative way, based 
on a simple visual inspection. The simulations make use 
of the value of r listed in Table 1, and vary u within its 
empirically defined limits for each quasar (Table 2). The 
flare energy and background luminosity are determined 
self consistently through eqs. (9) and (10). Square shaped 
flares were used for consistency with the results already 
presented. 

In Fig. 3 we present three illustrative examples of this 
experiment. The top curves in each panel show the B 
band light curves for PG 0052, 1309 and 1404, as observed 
in the Wise Observatory campaign. The thin solid line 
in the bottom illustrates a randomly chosen realization of 
the light curve, computed with v set to half-way between 
Vmin and Vmax- The thick solid line shows the simulation 
which, among thousands of runs, best matches the ob- 
served light curve in a sense, whereas the dashed curve 
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Fig. 2. — Relation between the net variability indices computed from the light curve, 5 = a/L, and those computed with the asymptotic 
variance inferred from the SF. (a) B-band. (b) R band, (c) Comparison of B and R band net variability. 



shows the same simulation after sampling it with the same 
rest-frame pattern as the data and with Gaussian pertur- 
bations added mimicking the sequence of photometric er- 
rors reported in G99. All model light curves were vertically 
displaced for clarity. By construction, their mean luminos- 
ity is essentially identical to the observed one, as are their 
variances. 

Considering the simplicity of the model and the require- 
ment of consistency with the constraints defined in §3.3, 
the resemblance between observations and the models is 
very satisfactory. The "best-matches" are particularly 
striking. Even better matches would be obtained fitting 
the light curves, i.e., adjusting both the global parameters 
and the individual flare birth-dates to optimize the residu- 
als. Of course, the meaning and usefulness of such detailed 
fits would be questionable given the number of parameters 
involved. For instance, there are anywhere between 36 and 
145 flares for the length of the PG0052 light curve, and up 
to 1316 in PG1404. Indeed, with so many degrees of free- 
dom, it might be impossible to ever disprove such models. 
Unlike the flare dates, however, one has much less liberty 
to play with the global parameters, and our estimates of 
T, v, E\ and c\ go a long way towards constraining the 
simulations to a region of parameter space capable of pro- 



ducing quasar looking light curves. Simulations with much 
longer or shorter time scales, for instance, do not resemble 
the observations at all. We therefore conclude that even a 
simple Poissonian superposition of square flares is capable 
of producing light curves which are very similar to quasars 
and simultaneously compatible with their basic light curve 
statistics. 

3.4.1. Structure Function 

The SFs of PG 0052 and 1404 are compared to the cor- 
responding simulated light curves in Fig. 4. As antici- 
pated, relatively large deviations from the smooth theoret- 
ical curve occur even for model light curves. The bottom 
panels show how the SF is improved with light curves twice 
to five times longer than the observed ones. These were 
constructed simply extending the sampling of the model 
by patching together 2 (panels e and f) or 5 (g and h) se- 
quences of the actual observing dates for these two quasars, 
thus maintaining the G99 sampling pattern. Even for such 
long (13-45 yr) hypothetical light curves the SF oscillates, 
though one sees it gradually converging to its statistically 
expected shape and amplitude. Increasing the sampling 
rate is not as benefic as increasing the length of the data 
train, as SF oscillations on intervals Ai are only averaged 
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Fig. 3. — Comparison of observed quasar B-band light curves (top dotted lines) and simulations. Thin solid lines show a randomly chosen 
simulation. The thick lines indicate a "best match" model (see text), while the dotted line show the same model after sampling as in the data. 
All light curves are normalized to Lb (Table 1), and the simulations are shifted downwards for clarity. The global Poissonian parameters for 
the simulations in thick hues are t = 0.70, 1.87 and 0.47 yr; u = 19, 47 and 223 yr"!; Eb = 4.8 X 10'''^, 1.8 X lO"'^ and 4.9 X 10*5 ergA"!; 
cg = 15, 8 and 7% for PG 0052, 1309 and 1404 respectively. The simulations in the bottom thin lines have the same lifetimes as above, but 



1/ = 16, 29 and 131 yr"!; Eb = 5.2 X 10*'^, 2.4 X 10*^ and 6.3 X lO"'^ ergA" 
These values are all, by construction, consistent with the limits in Table 2. 



cs = 21, 28 and 29% for PG 0052, 1309 and 1404 respectively. 



out for At 3> T, i.e, as independent (separated by ^ t) 
portions of the light curve are sampled a statistically sig- 
nificant number of times. This explains why even long 
light curves exhibit long term SF 'noise' but are much 
better defined on short time scales (Fig. 4e-h). 

These experiments illustrate the inherent difficulty in 
modeling SFs of stochastic processes. As a further exam- 
ple, we note that the SF analysis of a longer (11 yr) light 
curve of PG 1226 (= 3C 273) by Paltani, Courvoisier & 



Walter (1988) suggests a value of r of ~ 0.6 yr, a factor 
of 4 smaller than the 2.3 yr we found for the 5.5 yr G99 
monitoring of this quasar. The SF of PG 1226 is similar 
to that of PG 1229 in the B band (Fig. 1), with an initial 
peak at small At (~ 0.4 yr), followed by oscillations on 
longer time scales. Whereas for PG 1229 the global fits 
favored the small r, in PG 1226 the amplitude of the 0.5 
yr peak is smaller, and a longer r was favored, similar to 
what is seen in the R band SF of PG 0838 (Fig. 1). Even 
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Fig. 4. — Structure Function of PG 0052 (panel a) and PG 1404 (b) in the B band. Panels c and d show the SFs of the corresponding 
simulated light curves (shown as dashed lines in Fig. 3). Simulations twice and five times longer than the length of the Wise Observatory 
monitoring of these quasars were used to compute the SFs in panels e and f (2x), and g and h (5x). In all plots the solid line marks the same 
combined B and R SF fit performed upon the observed light curves. All SFs are normalized to twice the observed tr^ . 



considering the differences in technicalities of SF analy- 
sis between different works, one is forced to conclude that 
the actual uncertainties in the SF parameters is likely to 
substantially exceed the formal fit-errors (about 5% of the 
parameter values). 

Overall, we conclude that SF analysis for individual ob- 
jects is just starting to become possible, despite the excel- 
lent quality of the G99 light curves. We are however con- 
fident that the fits performed in this paper are correct to 
within better than an order of magnitude, and thoroughly 
meet our goal of providing rough but useful estimates for 



the basic parameters of Poissonian models for the optical 
variability of quasars. 

4. DISCUSSION 

Poissonian models provide a simple and elegant mathe- 
matical framework which relates basic parameters to mea- 
surable quantities, as demonstrated by the results pre- 
sented above. The same technique can be applied to other 
data sets and spectral bands, thus increasing the number 
of constraints. In particular, it will be interesting to apply 
it to the spectrophotometric monitoring observations of 
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28 quasars, all in the G99 sample, reported in Kaspi etal. 
(2000). This should yield a clearer picture of the (optical) 
spectral energy distributions of the variable and constant 
components, as well as better constrained parameters than 
was possible with only two wavelengths. An illustration of 
how important it is to properly define Lx and 6x is given 
in Fig. 5. 

The Eb limits listed in Table 2 encompass both the es- 
timates of Cid Fernandes etal. (1996) and Garcia etal. 
(1999) for different samples. Their estimates, however, 
were based on the properties of the ensemble of quasars, 
and assumed typical variability time-scales, whereas here 
we attempted to treat each quasar individually. In fact, 
this is the first time that constraints on Poissonian param- 
eters are computed on an object by object basis for a large 
number of quasars. 

A further comparison can be made with the results of 
Paltani etal. (1988) for 3C273. Wc estimate from their 
plots a monochromatic fiare energy of ~ 2 X ergA-\ 
about an order of magnitude smaller than EB,min in Table 
2. This discrepancy is explained by the already mentioned 
~ 5 times smaller r found by those authors, as well as 
their ^ 2 times smaller asymptotic net variability com- 
pared to the 6sF,B = 0.13 in Table 2, and another factor 
of 1.5 due to the redshift effect (see below). The dis- 
crepancy in the asymptotic SF variance is mostly due to 
the fact that they allow for a slowly (> 10 yr) varying 
component, attributed to a blazar behavior. Two of the 
other 6 radio loud objects in G99 sample (PG 1100 and 
1512) show evidence for long time-scale variations which 
could be due to a similar slow component (marginal evi- 
dence that the radio loud objects in this sample present 
longer time-scales can be seen in Fig. 6c). It is in fact 
adequate to recall that AGN variability has often been de- 
scribed in terms of the superposition of rapid and slowly 
varying components (e.g., Lyutyi 1977, Pica etal. 1988). 
A caveat in the present analysis is thus that what we have 
been calling a "constant component" may well correspond 
to an underlying slow process. 

We note that no attempt was made to correct any of 
the variability indices nor the inferred parameters for the 
variability-A effect which plague fixed band photometric 
light curves of quasars, leading to overestimated variabil- 
ity amplitudes (Cristiani etal. 1996; Cid Fernandes etal. 
1996; Aretxaga etal. 1997; Garcia etal. 1999). Fortu- 
nately, the low redshifts of the G99 quasars minimize this 
effect. For the median z of the sample, and using the pa- 
rameterization of the 6 X X relation of Garcia et al. (1999), 
wc estimate that the variability indices in Table 1 (columns 
4 to 7) would need to be multiplied by typically 0.8 to ob- 
tain the rest-frame indices. This would increase Mmin, 
■N'max and the corresponding u limits by some 56% while 
reducing the energies by a factor of 0.64. These are rela- 
tively small factors considering the allowed ranges for the 
parameters and were not applied in Table 2. Further- 
more, such corrections would soon become obsolete since 
the spectral data collected by Kaspi et al. (2000) allows the 
direct computation of rest-frame continuum variability in- 
dices. We did, however, experiment with the A-correction 
in the correlation analysis presented next. 

4.1. Correlation Analysis 



Perhaps the most puzzling result of G99 was the finding 
that there are few convincing correlations between vari- 
ability indices and a large array of other observed proper- 
ties for the PG quasars. In Table 3 we synthesize the re- 
sults of a correlation analysis analogous to that performed 
by G99, but for our SF-based variability indices. The 
table is similar to their Table 4, from which the multi- 
wavelength data was borrowed, and lists the percentage 
probabilities {Pr) of no correlation in a Spearman's rank 
test, small values indicating significant correlations and 
negative values corresponding to anticor relations. Entries 
between parentheses correspond to correlations after cor- 
rection of the variability amplitudes with the 5-\-z relation 
following Garcia etal. (1999). 

The anti-correlation between the variability amplitude 
(here measured by 5sf,b and 8sF,Fd and the optical lu- 
minosity found in other studies (e.g., Hook etal. 1994) is 
essentially absent for this sample, as illustrated in Fig. 6a. 
Somewhat more significant anti-correlations are obtained 
by applying the A correction (Fig. 6b; Table 3), but in 
both cases the correlations here are even weaker than those 
found by G99. The dotted lines in Fig. 6a and b mark the 
prediction of simple Poissonian models in which c\, E\ and 
T are identical for all quasars. This illustrates the often 
claimed failure of this model (e.g. Hook etal. 1994), and 
the necessity to allow for a diversity of values for these pa- 
rameters to account for the observed variability properties 
in the framework of Poissonian processes. 

The variability time-scale seems to increase with lumi- 
nosity (Fig. 6c), but, as in G99, with a large scatter. Previ- 
ous evidence for such a correlation was reported by Cris- 
tiani etal. (1996), in an analysis of the ensemble SF of 
486 objects spanning a much larger range of luminosities 
and redshifts. It is interesting to note that applying their 
Model E SF-fits to the typical absolute magnitudes of G99 
quasars yield r between ^ 0.5 and 2.5 yr, compatible with 
our individual SF fits. We also find r to correlate posi- 
tively with the 5sf's, but this may be due to a trade-off 
effect in the SF fitting (§3.4.1). The weak anticorrelation 
with Uox found by G99 is much stronger here (Fig. 6d). 

The good correlations between the variability amplitude 
and the equivalent widths of H/3 and Hell A 4686 found by 
G99 are reproduced at about the same high level of signif- 
icance (Fig. 7a). Unlike G99, however, we find the optical 
to X-ray index 

Oiox f*^ significantly anticorrclated with 
both 6sF,B (Fig. 7d) and 5sf,r- All these correlations are 
improved with the A-correction. 

A correlation was found between the rati o of the B to 
R asymptotic standard deviations and Lb as well as z. 
This could be due to the redshift effects already discussed, 
but in this case a 6-X law different from that of Garcia 
etal. (1999) would be implied, since the ^-effects upon 
cb,sf/o'k,sf cancels out in their parameterization. Alter- 
natively, if this interpretation is proved wrong, one would 
conclude that the flares are bluer in more luminous sources 
(eq. 9). 

Overall, we largely confirmed G99 results, with a few 
minor differences (e.g.. Fig. 6d) undoubtedly due to the 
different methodologies employed in the definition of vari- 
ability amplitudes and time-scales. Despite the few good 
correlations found, whose meaning we discuss below, a 
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Fig. 5. — Background components {Cx, dotted lines) for a hypothetical quasar with Lx oc A ^-^ (top solid line) and <5>, oc A ^. The plot 
illustrates how the spectral energy distribution of the constant component can be deduced from (observationally) well determined Lx and Sx, 
coupled with variability derived constraints upon cx at any particular wavelength. 



striking result of this analysis is the large scatter present 
even in the best correlations obtained. In any scenario 
for AGN, the variability properties must somehow be dic- 
tated by the physical conditions prevalent in the active 
nucleus (accretion rate, black hole mass, orientation, etc.). 
Since these conditions must also define other aspects of 
AGN phenomenology, it is natural to expect that proper- 
ties driven by common causes should exhibit some degree 
of correlation. Given the controversial history of correla- 
tions in quasar variability studies, however, it is perhaps 
not so surprising that so few good correlation were found. 

Regarding Poissonian models, one would obviously like 
to search for correlations between the basic parameters 
and other properties. The constraints upon A/", v, E\ and 
c\ derived in this paper are however not strong enough 
to warrant a proper correlation analysis. Taking the mid- 
dle value between Mmin and Afmax as a measure of A/", 
we obtain an anticorrelation between M and EW(H/3) and 
a positive correlation with aox^ both with at the 1% 
level. The flare rate follows the same trends. For the rea- 
son discussed above, it would be premature to give much 
emphasis to these correlations. 

4.1.1. Interpretation of the correlations with EW(HP) 

A simple interpretation of the significant correlation be- 
tween EW(H/3) and the variability amplitude is possible 
by postulating that the variable component dominates the 
ionization of the gas. This idea is also compatible with our 
earlier conclusion (§2) that the variability-wavelength an- 
ticorrelation indicates the dilution of a blue variable com- 
ponent by a red underlying background. In this hypothe- 
sis both variable (Vx) and constant {C\) components con- 



tribute to the continuum under H/?, but only the former 
is proportional to the ionizing luminosity (Lion), so one 
predicts 



EW{H(3) cx 



i^Eb 



vEf 



Ce 



(13) 



where we used the B band because of its proximity to II/3. 
This relation may be rewritten as 



EW{H(3) cx Af^/'^S^ 



SF,B 



Vb = I - Cb 



(14) 



which reveals a proportionality between EW(II/3) and 
^SF,B- Obviously, the same prediction applies to 
EW(Hell). 

Using TV = {J\fmin +J^niax)/2, the product M^/^SsF,B 
correlates at the Pr — 0.6% level with EW(II/3), in agree- 
ment with the prediction. Since our estimates of A/" are not 
independent of Ssf.b (eq. 11), it is perhaps more mean- 
ingful to test the predicted relation using cb, for which 
we have a robust upper limit imposed by the minimum in 
the light curve (/is in Table 1). In Fig. 7b we see that 
the expected anticorrelation is confirmed with a high sig- 
nificance (Table 3) . Even considering that fiB is an upper 
limit to and the scatter in the plot, it is interesting to 
see that the trend is roughly linear, as predicted. 

^B is also strongly correlated with aox (Pr — 0.1%), 
the X-ray spectral index ax {Pr = 1-8%), and anticorre- 
lated with EW([OIII]) {Pr = 0.8%), [OIII] to H/3 peak 
intensities ratio {Pr = 0.2%) and EW(Hell) {Pr = 0.5%). 
The correlation with aox, in particular, gives strength to 
our working hypothesis, insofar as this index can be inter- 
preted as indicative of the ratio between the constant and 
variable (ionizing) components. 
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Fig. 6. — Correlations between variability properties and luminosity (a— c) and aox (d). Filled circles mark radio quiet objects, whereas 
squares correspond to radio loud sources Dotted lines in the top panels are simple Poissonian models with 0.51og(l — cg)Eg/T = 40.5 (top 
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as obtained from the SF fits. 



By admitting that different quasars have different back- 
ground fractions, and that Vx dominates the ionization, 
we naturally explain at least some of the correlations ob- 
served. This argues against C\ being part of the ionizing 
source, such as hot but non- flaring portions of the accre- 
tion disk surface. Outer, colder regions of the disk can 
contribute to C\ , as do the host galaxy starlight and other 
circum-nuclear sources. 

The alternative hypothesis that Lion is proportional to 
the total V\ + C\ components would not induce the corre- 
lations above, while the hypothesis that Lion is governed 



by the constant component alone can be straightforwardly 
rejected, since this would not produce emission line vari- 
ability. 

We note that our eq. (13) is formally and philosophically 
identical to eq. (10) of Aretxaga & Terlevich (1994), who 
identify Cx with a young nuclear star-cluster and V\ with 
compact SN remnants exploding out of the same cluster 
(see also Aretxaga etal. 1997). Their analysis of EW(H/3) 
aimed to understand the near universal value of this quan- 
tity across the AGN luminosity scale (Binette, Fosbury & 
Parker 1993). In this model, a range of EW(H/3) values 



Quasar variability and Poissonian models 



13 



could be linked to either different SN explosion energies, 
or, as is more likely, to different H/3 and continuum pro- 
duction efficiencies in the SN remnants (see Cid Fernan- 
des 1997 for a review of the pros and cons of the starburst 
model) . Stochastic effects could also play a role in defining 
the EW(H/3) and variability relation, something which can 
happen in any Poissonian scenario. Indeed, the tentatively 
identified anticorrelation between M (and v) and EW(H/3) 
may be indicative of stochastic effects in the G99 data, and 
may be partly responsible for the scatter in Fig. 7a and d. 

4.2. Relation to physical models 

The ultimate goal of the Poissonian description of 
quasar variability is of course to extract information on 
the fundamental properties of the variability phenomenon 
and use them to guide and discriminate between physical 
theories. A detailed discussion of the implications of the 
estimates presented in §3 to different scenarios for AGN 
variability is beyond the scope of the present paper, but 
we would like to take one particular model as an illustra- 
tion of the possible associations between the observation- 
ally constrained parameters discussed here and physically 
meaningful quantities. 

In the model proposed by Haardt etal. (1994), blobs 
emerge from the surface of an accretion disk and release 
magnetically stored energy in the form of rapid X-ray 
flares. In the X-ray regime, the background fraction in 
this model would be associated with the fraction of the 
disk area covered by blobs (sec Figure 1 of Galeev, Rosner 
& Vaiana 1979), which equals the ratio of charge and dis- 
charge times. As they postulate the optical-UV variability 
to be driven by the X-ray flares, c\ is related to the fraction 
of the reprocessed luminosity and the optical-UV emission 
from the underlying stable component of the disk, but note 
that comparison with the limits on c\ derived in this pa- 
per require allowance for the diluting effects of larger scale 
sources, such as the host galaxy or an extended scatter- 
ing region. One of the predictions of this model is that at 
any time there are of order 10 active blobs, independent 
of the source luminosity, consistent with the limits on A/" 
established here. Also the reprocessed luminosity, which 
we would identify with {1 — c)L = vE in our formalism, 
is predicted to be similar to the X-ray luminosity. It is 
not clear, however, if the time smearing of the minutes 
long flares during the disk illumination which results in 
the reprocessing of the blobs energies can be made com- 
patible with the ~ 0.5-3 yr time scales implied by the 
SF analysis. Surely this must impose some limits upon 
the geometry of the reprocessing surface. Furthermore, it 
remains to be established whether the flare energies can 
be made compatible with the limits laid out in §3.3. A 
more detailed scrutiny of the Haardt etal. (1994) model 
must therefore await more theoretical developments. On 
the observational side, a Poissonian analysis of X-ray light 
curves of quasars (similar to that done by Cid Fcrnandes 
1995 for EXOSAT light curves of Seyferts) could also set 
constraints upon this particular model. 

Despite their generality and appealing aspects, Poisso- 
nian models are by no means the only possible description 
of AGN variability. For instance, Mineshige & Shields 
(1990) showed that thermal limit cycles similar to those 
known to occur in dwarf novae (Warner 1995) can also 



take place in AGN disks, leading to eruptions followed 
by quiescent periods. Such a model, in which the degree 
of variability is "self-regulated" , would invalidate the for- 
malism employed here, which fundamentally rests upon 
the hypothesis of independence of the flares and statistical 
stationarity. The same applies to scenarios involving ran- 
dom walk-like or other kinds of state dependent behavior 
(e.g., Begelman & De Kool 1991; Stern, Svensson & Sikora 
1991) or periodic modulations associated to, for instance, 
precessing jets (e.g.. Abraham & Romero 1999). Test- 
ing the Poissonity of AGN light curves is however a hard 
task. Fourier phase coherence studies like that performed 
by Krolik, Done & Madejski (1993) for X-ray light curves 
of Seyferts, or intermittency tests (Vio et al. 1992) like the 
one carried out by Longo etal. (1996) for the historical 
light curve of NGC 4151, can in principle help discrimi- 
nating between Poissonian models and scenarios where the 
variations are due to coherent oscillations of a single en- 
tity, but these require more data points than are currently 
available from quasar optical monitoring studies. There- 
fore, until proven wrong, the Poissonian description may 
be regarded as a useful tool to unveil physically meaningful 
properties of AGN. 

5. SUMMARY AND CONCLUSIONS 

We have reviewed the Poissonian formalism for quasar 
variability in an attempt to provide a general framework 
which allows fundamental parameters to be estimated 
from good monitoring data. This was applied to the 6 
yr long B & R light curves of 42 PG quasars obtained 
by G99, yielding constraints for the energy, rate and life- 
times of the flares. Our main results can be summarized 
as follows. 

(1) The only reasonable way to account for the fact that 
quasars vary more at shorter wavelengths within a Poisso- 
nian scenario is to include the diluting effects of an under- 
lying "non-variable background" component redder than 
the spectral energy distribution of the putative flares. 

(2) A wide range of flare energies, lifetimes, and/or back- 
ground fractions has to be invoked to account for the ob- 
served variability properties (amplitudes and time-scales). 
This "stretching" of the simplest Poissonian scenario (in 
which all parameters are the same for every quasar) is 
warranted by the model independent result that quasars 
present a wide range of variability time scales, from ^ 0.5 
to more than 3 yr, as inferred from a Structure Function 
analysis. 

(3) The mean number of living flares is constrained to 
be of order M ~ 5 to 100, and lower limits for the R- 

band background contribution of typically 25% are estab- 
lished. These estimates are independent of cosmology, K- 
correction and extinction, and little sensitive to the SF 
analysis. 

(4) Flare rates between ~ 1-100 yr~^ and monochro- 
matic flare energies in the ^ ]^q46-48 Qj-g^-i range are 
implied by the data. Overall, the Poissonian parameters 
for individual quasars are constrained to within about an 
order of magnitude. 

(5) Light Curve simulations were performed and demon- 
strate the ability to reproduce the observed morphology 
of quasars light curves extremely well even for a blatantly 
simplistic "on/off" square-shaped model for the evolution 
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Fig. 7. — Correlations between Ssf,B) IJ-b, EW(H/3), EW(Hell) and aox (the larger oiox the softer the spectrum). Note that fiB, the fraction 
of the minimum to the mean B -band luminosity, is to be seen as an upper limit to true fraction of the background luminosity {cb < /^b)- 
Symbols as in Fig. 6. 



of the individual flares. Indeed, the whole Poissonian anal- 
ysis is highly insensitive to the flare shape. 

(6) Experiments were performed on the use of higher 
moments of the light curve as further constraints, but 
found to be of little use at present. 

(7) The variability properties of the PG quasars present 
little correlation with other properties. Even the best cor- 
relations identified, like that between the variability am- 
plitude and EW(H/3), present a substantial scatter, con- 
firming the results of G99. 

(8) The EW(H/3) x variability amplitude was inter- 
preted in a scenario where only the variable component 
participates in the ionization of the line emitting gas, 
consistent with conclusion (1) above. Correlations with 
EW(Hell) and the X-ray to optical spectral index further 
support this interpretation. 

Progress on this line of approach to AGN variability will 



require longer light curves to obtain more accurate esti- 
mates of the variability properties, which are important to 
constrain both Poissonian and non-Poissonian models. Re- 
gardless of possible improvements on the time-series anal- 
ysis techniques, it is important to encourage the contin- 
uation of the current CCD monitoring campaigns. Spec- 
tral information, some of which is already available (Kaspi 
etal. 2000), will also be valuable in deriving better con- 
straints for the model parameters and a more detailed pic- 
ture of the spectral behavior of the flares and background 
component. 

The most appealing aspect of Poissonian models is their 
generality. Scenarios as diverse as accretion disk instabil- 
ities, stellar collisions and supernovae all fall under the 
large Poissonian "umbrella". Detailed modeling will be 
required to bridge the gap between this mathematical for- 
malism and the physics of AGN variability and to fully 
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explore the constraints made possible by the application 
of this technique. 
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APPENDIX 

SKEWNESS AND KURTOSIS IN POISSONIAN MODELS 

The Poissonian formalism allows the computation of higher moments of the L{t) distribution, though most applications 
to AGN variability stop on the variance (eq. [7]). In this appendix we present expressions for the theoretically expected 
third (skewness) and fourth (kurtosis) moments, and discuss their applicability to the G99 data set. These moments were 
not used to constrain the flare properties of quasars in the main text, and are presented here only for the purposes of 
completeness and future reference. 

The skewness (7) and kurtosis (k) are defined by 

7 = a-3[L-L]3 (Al) 



K = (7-^[L - LY - ^ (A2) 

The —3 term in (A2) makes the kurtosis of a Gaussian distribution 0. Expressions for 7 and k can be derived extending 
the formalism employed by Cid Fernandes (1995), which yields (Vieira 2000) 
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7 = M^r)-'/' (A3) 

« = Mi't)-^ (A4) 
In these expression and play the role of "shape factors" : 

h = ^J I'm (A5) 



(A6) 



Both (py and equal 1 for square shape flares. Exponential flares result in (/)^ = 4/3 and = 2, while for triangular flares 

(py = 9/8 and (p^ = 27/20. The above relations reveal interconnections between the different moments in a Poissonian 
model: 7 oc (5 and k (x 6^, where S — a/L (eq. [7]). Unlike for d, an underlying non- variable component docs not affect 7 
nor K. In principle, this allows an estimate of the background fraction from the ratio of 7 (or k) to 5. 

The main problem concerning the use of higher light curve moments is that, as can be seen in eqs. (A3) and (A4), the 
predicted moments rapidly tend to the Gaussian limit (7 = k = 0) with increasing J\f = vt, i.e., as the superposition of 
events grows higher. This is obviously a consequence of the Central Limit Theorem (e.g., Papoulis 1965). 

The detection of deviations from Gaussianity in quasar light curves is known to be problematic (Press & Rybicki 1997), 
and the G99 light curves are no exception. In Fig. A8 we present the skewness and kurtosis for the Wise Observatory 
data. Corrections for the effects of photometric errors upon the moments were applied, but were negligible, given the 
excellent acciiracy of the G99 photometry. At first sight, the results in Fig. A8 would seem to immediately rule out any 
Poissonian model, as one finds negative 7 and k for about half of the objects, whilst the theory (eqs. [A3] and [A4]) predict 
only positive values! 

This, however, is likely to be an effect of sampling. To demonstrate this, we have run a series of Poissonian light curve 

simulations and compared the output 2nd, 3rd and 4th moments with the predicted values as a function of observational 
parameters such as the length of the light curve {Tobs) and the number of observations (riots). The results clearly show 
that negative 7 and k occur very often also in simulated light curves, particularly for Tobs % lOr, as illustrated in Fig. A9. 
The scatter decreases for increasing riobsi but the agreement between analytical and simulated moments is only achieved 
for large riobs and Tobs ^ t. This bias is essentially insensitive to flare proflle or their rate. The conclusion here is that 
one cannot use the currently available data to strongly constrain the higher moments of the L{t) process. On the positive 
side, these experiments reinforce our conclusion that the actual flare shape is irrelevant for most of the analysis presented 
in this paper (see also next section). 

We flnalize by noting that, in analogy with what was done for the second moment, higher moments should be computed 
from their asymptotic {At 00) behavior estimated via SFs of the corresponding order. These, however, are subjected 
to large uncertainties due to the high powers involved and are not presented here. 

STRUCTURE FUNCTIONS FOR SIMPLE FLARE PROFILES 
The structure function of a Poissonian sequence of flares is: 

SF{At) = SF{oo) s{At) (Bl) 

where SF{oo) = 2uE'^/t is twice the asymptotic variance of the L{t) process, and s{At), the normalized SF of individual 
flares, is given by 

and is sensitive only to the shape of flare time profile l{t). 

Square flares 

For l{t) = lo between Q <t <T and otherwise one obtains a a linear SF up to At = T, with 

siAt) = ^ (B3) 

and 1 for At > T. The life-time t (given by eq. 5) of square flares equals T. The dotted lines in Fig. BIO show the 
resulting SF. 
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Fig. A8. — Skewness and kurtosis for the G99 B-band light curves plotted against each other (top panel), and Sgp (bottom panels). Error 
bars for 7 and ft were computed via bootstrap, resampling riQij^ randomly chosen points (with repetition) from the light curve 1000 times and 
computing the standard deviation of the resulting moments. The solid lines show the theoretical predictions for 0-y = 0k = 1 (square flares). 
In the bottom panels the two lines correspond predictions for a = 50% background fraction (upper curves) and no background (c^ = 0, 
bottom curves). Note that c affects only the abscissa in these plots, while the K-7 relation is independent of c. 



Exponentially decaying flares 

For l{t) ~ Iqc^*-/'^ flares one finds 

s(Ai) = 1 - e-'^*/^ (B4) 
while the life-time is t = 2T. This SF is shown as a dot-dashed line in Fig. BIO. 

Symmetric Triangular flares 
Linear flares with equal rise and decay times 



Ut\ -I X / (1 + ; for -T < i < , . 



produce the following SF: 
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Fig. A9. — Fractional variability (5), skewness (7) and kurtosis (k) for simulations as a function of the length of the light curve in units of 
T. Dashed, dotted and solid lines correspond to riots = H) 101 and 1001 observations respectively. The middle curves indicate the median 
over 1000 light curve simulations, whereas the bottom and top lines mark the 16 and 84% percentiles respectively, such that 68% of the points 
lie between them. Square flares and a rate of 30 events per r were used in the simulations. The theoretical moments are indicated by the 
arrows on the right side of each panel. Large (and systematic in the cases of 5 and k) deviations from the expected values of the moments 
are observed for small n^i,^ and specially for short light curves {T^bs ^ lOr). 
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The relation between r and T for triangular shots is r = 3T/2. The solid line in Fig. BIO shows this SF. 

Asymmetric Triangular flares 
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Fig. bio. — Theoretical Structure Functions for flares with square exponential, triangular and asymmetric triangle profiles. 

Triangular flares with unequal rise (T^) and decay (T^) times 



produce more complicated SFs: 



(B7) 
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(B8) 



where T = Tr + Td. The sohition above corresponds to T,. < Td- The corresponding expression for flares that spend more 
time rising than decaying is identical, with T,. swapped by Td- The hfc-time in either case is given by t = 3r/4. G99 
finding that quasars spend more time on the rise tlian fading gives some motivation to using Tr > Td, but the SF, being 
a mean of squared differences, does distinguish between and (see Kawaguchi et al. 1998 for techniques to expfore 
the time-assymetry of flares). Fig. BIO shows the SF for a ratio of 10 between and (or vice-versa) as a dashed line. 



